# Load required libraries
library(rnaturalearth)
library(dplyr) # for the filter() function
library(readr)  # for read_csv()
library(RColorBrewer)
library(scales)
library(readxl)
suppressPackageStartupMessages(library(sp))
suppressPackageStartupMessages(library(sf))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(ggiraph))
suppressPackageStartupMessages(library(plotly))
library(htmlwidgets)

workDir <- "/Users/meijian/Work/SIMPLE20181102/"
setwd(workDir)

# Get world map data
africa_map <- ne_countries(continent = "Africa", returnclass = "sf")
africa_map$sovereignt[africa_map$sovereignt == "eSwatini"] <- "Eswatini"
africa_map$sovereignt[africa_map$sovereignt == "Ivory Coast"] <- "Côte d'Ivoire"

# grid info
grid_list <- paste0(workDir, "Input/Map/SIMPLE-All_Africa_Countries_List.xlsx")
grid_info <- read_excel(grid_list)
grid_info$Country[grid_info$Country == "Guinea Bissau"] <- "Guinea-Bissau"
grid_info$Country[grid_info$Country == "Swaziland"] <- "Eswatini"
grid_info$Country[grid_info$Country == "Ivory Coast"] <- "Côte d'Ivoire"
grid_info$Country[grid_info$Country == "S. Sudan"] <- "South Sudan"
country_list <- unique(grid_info$Country)
country_list <- country_list[!country_list %in% c("Egypt", "Djibouti")] # remove Egypt and Djibouti due to 100% irrigation

# capitalize first letter function
capitalize_first_letter <- function(string) {
  paste0(toupper(substr(string, 1, 1)), substr(string, 2, nchar(string)))
}

# yield_upper <- c()
yield_upper_fao <- c(1000, 20000, 2000, 1000, 2000, 4000, 12000, 1000, 2000, 10000, 1000, 800, 1000, 3000, 8000, 10000, 10000, 15000, 15000)
yield_upper_sim <- c(20000, 1000, 20000, 12000, 2000, 1000, 1000, 4000, 2000, 12000, 4000, 4000, 3000, 12000, 1000, 2000, 10000, 1000, 800, 1000, 3000, 8000, 10000, 1500, 10000, 15000)

production_upper_fao <- c(20, 40000, 2000, 80, 2000, 10000, 2000, 200, 6000, 25, 300, 4000, 600, 2000, 3000, 250, 1600, 30000)
production_upper_sim <- c(100, 20, 40000, 100, 2000, 150, 80, 300, 2000, 100, 1500, 10000, 40, 2000, 5000, 200, 6000, 25, 300, 4000, 600, 2000, 3000, 2000, 250, 30000)

#============================================================
# Obs plots
# Obs file Path
obs_file_path <- paste0(workDir, "Observation/FAO_Africa/yield/")

all_files <- list.files(path = obs_file_path, pattern = "\\.csv$")

africa_subset <- c()
species_list_FAO <- c()

for (i in 1:length(all_files)) {
  # Crop Name
  crop_name <- substr(all_files[i],1,nchar(all_files[i])-4)
  
  # Read the CSV file
  yield_obs <- read_csv(paste0(obs_file_path, all_files[i]),show_col_types = FALSE)
  
  # Summarize the value by country name, calculating the mean
  yield_obs_mean <- yield_obs %>%
    group_by(Area) %>%
    filter(Year >= 2000) %>%
    summarize(mean_value = mean(Value, na.rm = TRUE)/10)
  yield_obs_mean$Area[yield_obs_mean$Area == "Congo"] <- "Republic of the Congo"
  yield_obs_mean <- yield_obs_mean[!(yield_obs_mean$Area %in% c("Egypt", "Djibouti")), ] # remove Egypt and Djibouti due to 100% irrigation
  
  # convert tomato fresh weight to dry weight
  if (crop_name == 'tomato') {yield_obs_mean$mean_value <- yield_obs_mean$mean_value * 0.15}
  
  # Join the data
  africa <- left_join(africa_map, yield_obs_mean, by = c("sovereignt" = "Area"))
  colnames(africa)[colnames(africa) == 'mean_value'] <- "Yield"
  colnames(africa)[colnames(africa) == 'sovereignt'] <- "Country"
  africa$Yield <- round(africa$Yield, digits = 2)
  
  # Create the map
  p <- ggplot(data = africa) +
    geom_sf(aes(fill = Yield)) +
    # scale_fill_gradientn(name = "Yield\n(kg/ha)", colors = brewer.pal(9, "BuGn"), #limits = c(-100, 100),
    #                      oob = squish, guide = guide_colourbar(barwidth = 1.5, barheight = 10),
    #                      na.value = "transparent") +
    scale_fill_viridis_c(guide = guide_colorbar(title = "Yield (kg/ha)", barwidth = 1.5, barheight = 10), 
                         na.value = "transparent", limits = c(0, yield_upper_fao[i]), oob = oob_squish, direction = -1) +
    # theme_bw() + 
    theme_void() +
    theme(plot.title = element_text(size = 16, color = "blue", face = "bold", hjust = 0),
          legend.title = element_text(size = 14), legend.text = element_text(size = 14))
    #labs(title = capitalize_first_letter(crop_name)) +
    # theme(
    #   text = element_text(size = 20),  # base font size, will affect all text elements
    #   plot.title = element_text(size = 20, color = "black", face = "bold", hjust = 0),  # font for titles
    #   axis.title = element_text(size = 16),  # font size for axis titles
    #   axis.text = element_text(size = 16),  # font size for axis text
    #   legend.title = element_text(size = 16),  # font size for legend title
    #   legend.text = element_text(size = 16),  # font size for legend items
    #   panel.grid.major = element_blank(),
    #   panel.grid.minor = element_blank()
    # )
  
  # Save the plot
  ggsave(filename = paste0("./Simple_output_discover/26crops/country_figures/yield/", crop_name, "_obs.png"), plot = p, width = 7, height = 6, units = "in", dpi = 300)
  
  ## interactive plot1
  africa_ia <- crosstalk::SharedData$new(africa)
  p <- ggplot(africa_ia) + geom_sf(aes(tooltip = Country, fill = Yield)) +
    # scale_fill_gradientn(name = "Yield\n(t/ha)", colors = brewer.pal(9, "BuGn"), #limits = c(-100, 100),
    #                      oob = squish, guide = guide_colourbar(barwidth = 1.5, barheight = 10),
    #                      na.value = "transparent") +
    scale_fill_viridis_c(guide = guide_colorbar(title = "Yield (t/ha)", barwidth = 1.5,
                                                barheight = 10), na.value = "transparent", direction = -1) +
    theme_void() +
    theme(plot.title = element_text(size = 16, color = "blue", face = "bold", hjust = 0),
          legend.title = element_text(size = 14), legend.text = element_text(size = 14))
    # labs(title = paste0(capitalize_first_letter(crop_name), " (FAO)")) +
    # theme(
    #   text = element_text(size = 20),  # base font size, will affect all text elements
    #   plot.title = element_text(size = 20, color = "black", face = "bold", hjust = 0.5),  # font for titles
    #   axis.title = element_text(size = 16),  # font size for axis titles
    #   axis.text = element_text(size = 16),  # font size for axis text
    #   legend.title = element_text(size = 16),  # font size for legend title
    #   legend.text = element_text(size = 16)  # font size for legend items
    # )
  
  plot <- ggplotly(p) %>%
    highlight(
      "plotly_hover",
      selected = attrs_selected(line = list(color = "black"))
    )
  
  #print(plot)
  # save interactive plot
  saveWidget(plot, file=paste0("./Simple_output_discover/26crops/country_figures/yield/", crop_name, "_obs.html"))
  
  africa_subset <- c(africa_subset, list(subset(africa$Country, !is.na(africa$Yield))))
  species_list_FAO <- c(species_list_FAO, crop_name)
}

#=====================================================================
## production plots
obs_file_path <- paste0(workDir, "Observation/FAO_Africa/production/")

all_files <- list.files(path = obs_file_path, pattern = "\\.csv$")
for (i in 1:length(all_files)) {
  # Crop Name
  crop_name <- substr(all_files[i],1,nchar(all_files[i])-4)
  
  # Read the CSV file
  yield_obs <- read_csv(paste0(obs_file_path, all_files[i]),show_col_types = FALSE)
  
  # Summarize the value by country name, calculating the mean
  yield_obs_mean <- yield_obs %>%
    filter(Value != 0) %>%
    filter(Year >= 2000) %>%
    group_by(Area) %>%
    summarize(mean_value = mean(Value, na.rm = TRUE))
  yield_obs_mean$Area[yield_obs_mean$Area == "Congo"] <- "Republic of the Congo"
  yield_obs_mean <- yield_obs_mean[!(yield_obs_mean$Area %in% c("Egypt", "Djibouti")), ] # remove Egypt and Djibouti due to 100% irrigation
  
  # convert tomato fresh weight to dry weight
  if (crop_name == 'tomato') {yield_obs_mean$mean_value <- yield_obs_mean$mean_value * 0.15}
  
  # Join the data
  africa <- left_join(africa_map, yield_obs_mean, by = c("sovereignt" = "Area"))
  colnames(africa)[colnames(africa) == 'mean_value'] <- "Yield"
  colnames(africa)[colnames(africa) == 'sovereignt'] <- "Country"
  africa$Yield <- round(africa$Yield/1000, digits = 2)
  
  # Create the map
  p <- ggplot(data = africa) +
    geom_sf(aes(fill = Yield)) +
    # scale_fill_gradientn(name = "Production\n(t)", colors = brewer.pal(9, "BuGn"), #limits = c(-100, 100),
    #                      oob = squish, guide = guide_colourbar(barwidth = 1.5, barheight = 10),
    #                      na.value = "transparent") +
    scale_fill_viridis_c(guide = guide_colorbar(title = "Production\n(kiloton)", barwidth = 1.5, barheight = 10), 
                         na.value = "transparent", limits = c(0, production_upper_fao[i]), oob = oob_squish, direction = -1) +
    # theme_bw() + 
    #theme_minimal() +
    theme_void()  +
    theme(plot.title = element_text(size = 16, color = "blue", face = "bold", hjust = 0),
          legend.title = element_text(size = 14), legend.text = element_text(size = 14))
    #labs(title = capitalize_first_letter(crop_name)) +
    # theme(
    #   text = element_text(size = 20),  # base font size, will affect all text elements
    #   plot.title = element_text(size = 20, color = "black", face = "bold", hjust = 0),  # font for titles
    #   axis.title = element_text(size = 16),  # font size for axis titles
    #   axis.text = element_text(size = 16),  # font size for axis text
    #   legend.title = element_text(size = 16),  # font size for legend title
    #   legend.text = element_text(size = 16),  # font size for legend items
    #   panel.grid.major = element_blank(),
    #   panel.grid.minor = element_blank()
    # )
  
  # Save the plot
  ggsave(filename = paste0("./Simple_output_discover/26crops/country_figures/production/", crop_name, "_obs_prod.png"), plot = p, width = 7, height = 6, units = "in", dpi = 300)
  
  ## interactive plot1
  africa_ia <- crosstalk::SharedData$new(africa)
  p <- ggplot(africa_ia) + geom_sf(aes(tooltip = Country, fill = Yield)) +
    scale_fill_gradientn(name = "Production\n(t)", colors = brewer.pal(9, "BuGn"), #limits = c(-100, 100),
                         oob = squish, guide = guide_colourbar(barwidth = 1.5, barheight = 10),
                         na.value = "transparent") +
    #theme_minimal() +
    theme_void() +
    labs(title = paste0(capitalize_first_letter(crop_name), " (FAO)"))  +
    theme(plot.title = element_text(size = 16, color = "blue", face = "bold", hjust = 0),
          legend.title = element_text(size = 14), legend.text = element_text(size = 14))
    # theme(
    #   text = element_text(size = 20),  # base font size, will affect all text elements
    #   plot.title = element_text(size = 20, color = "black", face = "bold", hjust = 0.5),  # font for titles
    #   axis.title = element_text(size = 16),  # font size for axis titles
    #   axis.text = element_text(size = 16),  # font size for axis text
    #   legend.title = element_text(size = 16),  # font size for legend title
    #   legend.text = element_text(size = 16)  # font size for legend items
    # )
  
  plot <- ggplotly(p) %>%
    highlight(
      "plotly_hover",
      selected = attrs_selected(line = list(color = "black"))
    )
  
  #print(plot)
  # save interactive plot
  saveWidget(plot, file=paste0("./Simple_output_discover/26crops/country_figures/production/", crop_name, "_obs_prod.html"))
}

#=====================================================================
## harvest area plots
obs_file_path <- paste0(workDir, "Observation/FAO_Africa/area/")

all_files <- list.files(path = obs_file_path, pattern = "\\.csv$")
for (i in 1:length(all_files)) {
  # Crop Name
  crop_name <- substr(all_files[i],1,nchar(all_files[i])-4)
  
  # Read the CSV file
  yield_obs <- read_csv(paste0(obs_file_path, all_files[i]),show_col_types = FALSE)
  
  # Summarize the value by country name, calculating the mean
  yield_obs_mean <- yield_obs %>%
    filter(Value != 0) %>%
    filter(Year >= 2000) %>%
    group_by(Area) %>%
    summarize(mean_value = mean(Value, na.rm = TRUE))
  yield_obs_mean$Area[yield_obs_mean$Area == "Congo"] <- "Republic of the Congo"
  yield_obs_mean <- yield_obs_mean[!(yield_obs_mean$Area %in% c("Egypt", "Djibouti")), ] # remove Egypt and Djibouti due to 100% irrigation
  
  # Join the data
  africa <- left_join(africa_map, yield_obs_mean, by = c("sovereignt" = "Area"))
  colnames(africa)[colnames(africa) == 'mean_value'] <- "Yield"
  colnames(africa)[colnames(africa) == 'sovereignt'] <- "Country"
  africa$Yield <- round(africa$Yield, digits = 2)
  
  # Create the map
  p <- ggplot(data = africa) +
    geom_sf(aes(fill = Yield)) +
    # scale_fill_gradientn(name = "Production\n(t)", colors = brewer.pal(9, "BuGn"), #limits = c(-100, 100),
    #                      oob = squish, guide = guide_colourbar(barwidth = 1.5, barheight = 10),
    #                      na.value = "transparent") +
    scale_fill_viridis_c(guide = guide_colorbar(title = "Harvested Area\n(ha)", barwidth = 1.5,
                                                barheight = 10), na.value = "transparent", direction = -1) +
    # theme_bw() + 
    #theme_minimal() +
    theme_void()  +
    theme(plot.title = element_text(size = 16, color = "blue", face = "bold", hjust = 0),
          legend.title = element_text(size = 14), legend.text = element_text(size = 14))
  #labs(title = capitalize_first_letter(crop_name)) +
  # theme(
  #   text = element_text(size = 20),  # base font size, will affect all text elements
  #   plot.title = element_text(size = 20, color = "black", face = "bold", hjust = 0),  # font for titles
  #   axis.title = element_text(size = 16),  # font size for axis titles
  #   axis.text = element_text(size = 16),  # font size for axis text
  #   legend.title = element_text(size = 16),  # font size for legend title
  #   legend.text = element_text(size = 16),  # font size for legend items
  #   panel.grid.major = element_blank(),
  #   panel.grid.minor = element_blank()
  # )
  
  # Save the plot
  ggsave(filename = paste0("./Simple_output_discover/26crops/country_figures/harvest_area/", crop_name, "_obs_area.png"), plot = p, width = 7, height = 6, units = "in", dpi = 300)
  
  ## interactive plot1
  africa_ia <- crosstalk::SharedData$new(africa)
  p <- ggplot(africa_ia) + geom_sf(aes(tooltip = Country, fill = Yield)) +
    scale_fill_gradientn(name = "Harvested Area\n(ha)", colors = brewer.pal(9, "BuGn"), #limits = c(-100, 100),
                         oob = squish, guide = guide_colourbar(barwidth = 1.5, barheight = 10),
                         na.value = "transparent") +
    #theme_minimal() +
    theme_void() +
    labs(title = paste0(capitalize_first_letter(crop_name), " (FAO)"))  +
    theme(plot.title = element_text(size = 16, color = "blue", face = "bold", hjust = 0),
          legend.title = element_text(size = 14), legend.text = element_text(size = 14))
  # theme(
  #   text = element_text(size = 20),  # base font size, will affect all text elements
  #   plot.title = element_text(size = 20, color = "black", face = "bold", hjust = 0.5),  # font for titles
  #   axis.title = element_text(size = 16),  # font size for axis titles
  #   axis.text = element_text(size = 16),  # font size for axis text
  #   legend.title = element_text(size = 16),  # font size for legend title
  #   legend.text = element_text(size = 16)  # font size for legend items
  # )
  
  plot <- ggplotly(p) %>%
    highlight(
      "plotly_hover",
      selected = attrs_selected(line = list(color = "black"))
    )
  
  #print(plot)
  # save interactive plot
  saveWidget(plot, file=paste0("./Simple_output_discover/26crops/country_figures/harvest_area/", crop_name, "_obs_area.html"))
}

#=====================================================================
# SIMPLE results, sim plots
species_list <- c("africaneggplant", "bambaragroundnut", "cassava", "cocoyam", "cowpea", "fingermillet", "fonio", 
                  "grasspea", "groundnut", "josephscoat", "lablab", "maize", "mungbean", "okra", "pearlmillet", "pigeonpea", 
                  "plantain", "pumpkin", "sesame", "sorghum", "soybean", "sweetpotato", "taro", "tef", "tomato", "yams")

speciesNum <- length(species_list)

filterByCROPGRIDS <- function(inPath, species, country_data) {
  library(ncdf4)
  library(raster)
  
  filtered_data <- data.frame()
  area_threshold <- 5
  mask_raster <- raster(xmn=min(country_data$long)-0.25, xmx=max(country_data$long)+0.25,
                        ymn=min(country_data$lat)-0.25, ymx=max(country_data$lat)+0.25,
                        res= 0.5, crs="+proj=longlat +datum=WGS84")
  
  HarvestAreaFile <- paste0(inPath, "Input/Map/CROPGRIDSv1.05_NC_maps/", "CROPGRIDSv1.05_", species, ".nc")
  
  if (file.exists(HarvestAreaFile)){
    nc_data <- nc_open(HarvestAreaFile)
    lon <- ncvar_get(nc_data, "lon")
    lat <- ncvar_get(nc_data, "lat")
    ha <- ncvar_get(nc_data, "harvarea")
    nc_close(nc_data)
    
    rownames(ha) <- lon
    colnames(ha) <- lat
    
    ha_df <- as.data.frame.table(ha)
    colnames(ha_df) <- c('lon','lat','ha')
    ha_df$lon <- as.numeric(as.character(ha_df$lon))
    ha_df$lat <- as.numeric(as.character(ha_df$lat))
    # convert Ocean = -2, Not assessed = -1 to 0
    ha_df$ha[ha_df$ha < 0] <- 0
    
    mask_ha <- rasterize(ha_df[,c('lon','lat')], mask_raster, ha_df[,'ha'], fun=sum, na.rm=T)
    filterTable_data <- data.frame(cbind(xyFromCell(mask_ha, 1:ncell(mask_ha)), round(values(mask_ha),2)))
    colnames(filterTable_data) <- c('long','lat','Area_sum')
    
    filter_data <- subset(filterTable_data, Area_sum >= area_threshold)
    
    filtered_data <- merge(country_data, filter_data,  by = c("long", "lat"), all.x = FALSE, sort=FALSE)
  }
  else {stop(paste0("AutoGenGridInFile.R: Can't find CROPGRIDS data for ",species,"!"))}
  
  return(filtered_data)
}

harvestArea <- list()
#yearRange <- c(1990:2019)
yearRange <- c(2000:2019) #baseline
#gcm <- c("GFDL", "IPSL", "MPI", "MRI")
gcm <- "ERA5" #baseline
out_summary <- data.frame(crop = as.character(NA), country = as.character(NA), gcm = as.character(NA), year = as.integer(NA), 
                          yield = as.numeric(NA), biomass = as.numeric(NA), duration = as.numeric(NA)) 
out_summary_prod <- data.frame(crop = as.character(NA), country = as.character(NA), gcm = as.character(NA), year = as.integer(NA), 
                          production = as.numeric(NA), biomass = as.numeric(NA), duration = as.numeric(NA)) 

for (species in species_list) {
  harvestArea <- filterByCROPGRIDS(workDir, species, grid_info)
  
  #pattern <- paste0(species, "_historical.*\\.csv$")
  pattern <- paste0('^', species, "_baseline.*\\.csv$") # make sure string start with species, no leading characters
  files <- list.files(path = paste0(workDir, "Simple_output_discover/26crops"), pattern = pattern, full.names = TRUE)
  sim_results <- list()
  for (file in files) {
    sim_results <- c(sim_results, list(read.table(file, header=TRUE, sep=",", stringsAsFactors=FALSE)))
  }
  #list.files(pattern = paste0(workDir, "Simple_output_discover/23crops/", species, "_historical*\\.csv$"))
  for (country in country_list) {
    harvestArea_country <- harvestArea[harvestArea$Country == country,]
    if (nrow(harvestArea_country) == 0) {next}
    else{
      for (i in 1:length(sim_results)) {
        sim_results_model <- sim_results[[i]]
        for (year in yearRange) {
          yield_ha_country <- merge(harvestArea_country, sim_results_model[sim_results_model$SowingYear == year, ], by = c("long", "lat"))
          out_onerow <- as.data.frame(t(c(species, country, gcm[i], year, sum(yield_ha_country$Area_sum * yield_ha_country$Yield) / sum(yield_ha_country$Area_sum), 
                                          sum(yield_ha_country$Area_sum * yield_ha_country$Biomass) / sum(yield_ha_country$Area_sum), 
                                          sum(yield_ha_country$Area_sum * yield_ha_country$Duration) / sum(yield_ha_country$Area_sum))))
          colnames(out_onerow) <- c("crop","country","gcm","year","yield","biomass","duration")
          out_summary <- rbind(out_summary, out_onerow)
          
          #production
          out_onerow_prod <- as.data.frame(t(c(species, country, gcm[i], year, sum(yield_ha_country$Area_sum * yield_ha_country$Yield), 
                                          sum(yield_ha_country$Area_sum * yield_ha_country$Biomass), 
                                          sum(yield_ha_country$Area_sum * yield_ha_country$Duration))))
          colnames(out_onerow_prod) <- c("crop","country","gcm","year","production","biomass","duration")
          out_summary_prod <- rbind(out_summary_prod, out_onerow_prod)
        }
      }
    }
  }
}
out_summary <- out_summary[-1,]
# write.csv(out_summary, file = "./Simple_output_discover/23crops_country.csv", row.names = FALSE)
# Read data from a CSV file into a dataframe
# out_summary <- read.csv(file = "./Simple_output_discover/23crops_country.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
out_summary$yield <- as.numeric(out_summary$yield)
out_summary$yield[is.nan(out_summary$yield)] <- NA
out_summary$biomass <- as.numeric(out_summary$biomass)
out_summary$biomass[is.nan(out_summary$biomass)] <- NA
out_summary$duration <- as.numeric(out_summary$duration)
out_summary$duration[is.nan(out_summary$duration)] <- NA
write.csv(out_summary, "./Simple_output_discover/26crops/country_figures/country_yield_base.csv", row.names=FALSE)

out_summary_prod <- out_summary_prod[-1,]
# write.csv(out_summary, file = "./Simple_output_discover/23crops_country.csv", row.names = FALSE)
# Read data from a CSV file into a dataframe
# out_summary <- read.csv(file = "./Simple_output_discover/23crops_country.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
out_summary_prod$production <- as.numeric(out_summary_prod$production)
out_summary_prod$production[is.nan(out_summary_prod$production)] <- NA
out_summary_prod$biomass <- as.numeric(out_summary_prod$biomass)
out_summary_prod$biomass[is.nan(out_summary_prod$biomass)] <- NA
out_summary_prod$duration <- as.numeric(out_summary_prod$duration)
out_summary_prod$duration[is.nan(out_summary_prod$duration)] <- NA

# plot maps
for (species in species_list) {
  sim_data_map <- data.frame(country = NA, yield = NA)
  out_summary_species <- out_summary[out_summary$crop == species,]
  for (country in unique(out_summary_species$country)) {
    out_sum_sp_ct <- out_summary_species[out_summary_species$country == country, ]
    out_sum_sp_ct_y <- mean(out_sum_sp_ct$yield, na.rm = TRUE)
    out_sum_sp_ct_y <- as.data.frame(t(c(country, out_sum_sp_ct_y)))
    colnames(out_sum_sp_ct_y) <- c("country", "yield")
    sim_data_map <- rbind(sim_data_map, out_sum_sp_ct_y)
  }
  sim_data_map <- sim_data_map[-1,]
  sim_data_map$yield <- as.numeric(sim_data_map$yield)
  
  # subset based on the countries with available FAO data
  tryCatch({
  sim_data_map_subset <- sim_data_map[sim_data_map$country %in% africa_subset[[which(species == species_list_FAO)]], ]
  }, error = function(error_message) {
    sim_data_map_subset <<- sim_data_map
  })
  
  # Join the data
  africa <- left_join(africa_map, sim_data_map_subset, by = c("sovereignt" = "country"))
  colnames(africa)[colnames(africa) == 'yield'] <- "Yield"
  colnames(africa)[colnames(africa) == 'sovereignt'] <- "Country"
  africa$Yield <- round(africa$Yield, digits = 2)
  
  # Create the map
  p <- ggplot(data = africa) +
    geom_sf(aes(fill = Yield)) +
    # scale_fill_gradientn(name = "Yield\n(t/ha)", colors = brewer.pal(9, "BuGn"), #limits = c(-100, 100),
    #                      oob = squish, guide = guide_colourbar(barwidth = 1.5, barheight = 10),
    #                      na.value = "transparent") +
    scale_fill_viridis_c(guide = guide_colorbar(title = "Yield (kg/ha)", barwidth = 1.5, barheight = 10), 
                         na.value = "transparent", limits = c(0, yield_upper_sim[which(species == species_list)]), oob = oob_squish, direction = -1) +
    # theme_bw() + 
    theme_void() +
    theme(plot.title = element_text(size = 16, color = "blue", face = "bold", hjust = 0),
          legend.title = element_text(size = 14), legend.text = element_text(size = 14))
    #labs(title = capitalize_first_letter(crop_name)) +
    # theme(
    #   text = element_text(size = 20),  # base font size, will affect all text elements
    #   plot.title = element_text(size = 20, color = "black", face = "bold", hjust = 0),  # font for titles
    #   axis.title = element_text(size = 16),  # font size for axis titles
    #   axis.text = element_text(size = 16),  # font size for axis text
    #   legend.title = element_text(size = 16),  # font size for legend title
    #   legend.text = element_text(size = 16),  # font size for legend items
    #   panel.grid.major = element_blank(),
    #   panel.grid.minor = element_blank()
    # )
  
  # Save the plot
  ggsave(filename = paste0("./Simple_output_discover/26crops/country_figures/yield/", species, "_sim.png"), plot = p, width = 7, height = 6, units = "in", dpi = 300)
}

# plot production maps
for (species in species_list) {
  sim_data_map <- data.frame(country = NA, production = NA)
  out_summary_species <- out_summary_prod[out_summary_prod$crop == species,]
  for (country in unique(out_summary_species$country)) {
    out_sum_sp_ct <- out_summary_species[out_summary_species$country == country, ]
    out_sum_sp_ct_y <- mean(out_sum_sp_ct$production, na.rm = TRUE)
    out_sum_sp_ct_y <- as.data.frame(t(c(country, out_sum_sp_ct_y)))
    colnames(out_sum_sp_ct_y) <- c("country", "production")
    sim_data_map <- rbind(sim_data_map, out_sum_sp_ct_y)
  }
  sim_data_map <- sim_data_map[-1,]
  sim_data_map$production <- as.numeric(sim_data_map$production) / 1000000
  
  # subset based on the countries with available FAO data
  tryCatch({
    sim_data_map_subset <- sim_data_map[sim_data_map$country %in% africa_subset[[which(species == species_list_FAO)]], ]
  }, error = function(error_message) {
    sim_data_map_subset <<- sim_data_map
  })
  
  # Join the data
  africa <- left_join(africa_map, sim_data_map_subset, by = c("sovereignt" = "country"))
  colnames(africa)[colnames(africa) == 'production'] <- "Production"
  colnames(africa)[colnames(africa) == 'sovereignt'] <- "Country"
  africa$Production <- round(africa$Production, digits = 2)
  
  # Create the map
  p <- ggplot(data = africa) +
    geom_sf(aes(fill = Production)) +
    # scale_fill_gradientn(name = "Production\n(t)", colors = brewer.pal(9, "BuGn"), #limits = c(-100, 100),
    #                      oob = squish, guide = guide_colourbar(barwidth = 1.5, barheight = 10),
    #                      na.value = "transparent") +
    scale_fill_viridis_c(guide = guide_colorbar(title = "Production\n(kiloton)", barwidth = 1.5, barheight = 10), 
                         na.value = "transparent", limits = c(0, production_upper_sim[which(species == species_list)]), oob = oob_squish, direction = -1) +
    # theme_bw() + 
    theme_void() +
    theme(plot.title = element_text(size = 16, color = "blue", face = "bold", hjust = 0),
          legend.title = element_text(size = 14), legend.text = element_text(size = 14))
    #labs(title = capitalize_first_letter(crop_name)) +
    # theme(
    #   text = element_text(size = 20),  # base font size, will affect all text elements
    #   plot.title = element_text(size = 20, color = "black", face = "bold", hjust = 0),  # font for titles
    #   axis.title = element_text(size = 16),  # font size for axis titles
    #   axis.text = element_text(size = 16),  # font size for axis text
    #   legend.title = element_text(size = 16),  # font size for legend title
    #   legend.text = element_text(size = 16),  # font size for legend items
    #   panel.grid.major = element_blank(),
    #   panel.grid.minor = element_blank()
    # )
  
  # Save the plot
  ggsave(filename = paste0("./Simple_output_discover/26crops/country_figures/production/", species, "_sim_prod.png"), plot = p, width = 7, height = 6, units = "in", dpi = 300)
}

#=======================================================================
# animated plots
library(plotly)
library(ggplot2)
library(gapminder)

species_list <- c("africaneggplant", "bambaragroundnut", "cassava", "cocoyam", "cowpea", "fingermillet", "fonio", 
                  "grasspea", "groundnut", "josephscoat", "lablab", "maize", "mungbean", "okra", "pigeonpea", "sesame",
                  "sorghum", "soybean", "sweetpotato", "taro", "tef", "tomato", "yams")
crop_type_list <- c("Vegetables", "Legumes", "Roots/Tubers", "Roots/Tubers", "Legumes", "Cereals", "Cereals", 
               "Legumes", "Oilseeds", "Vegetables", "Legumes", "Cereals", "Legumes", "Vegetables", "Legumes", "Oilseeds", 
               "Cereals", "Legumes", "Roots/Tubers", "Roots/Tubers", "Cereals", "Vegetables", "Roots/Tubers")

df_sum <- data.frame(Crop = as.character(NA), CropType = as.character(NA), Year = as.integer(NA), 
                     Yield = as.numeric(NA), Biomass = as.numeric(NA), Duration = as.numeric(NA)) 
for (i in 1:length(species_list)) {
  #sim_data_map <- data.frame(country = NA, yield = NA)
  species <- species_list[i]
  crop_type <- crop_type_list[i]
  out_summary_species <- out_summary[out_summary$crop == species,]
  out_sum_sp_year_mean <- list()
  for (year in yearRange) {
    out_sum_sp_year <- out_summary_species[out_summary_species$year == year,]
    out_sum_sp_year_mean <- c(out_sum_sp_year_mean, mean(out_sum_sp_year$yield, na.rm=TRUE), 
                              mean(out_sum_sp_year$biomass, na.rm=TRUE), mean(out_sum_sp_year$duration, na.rm=TRUE))
  }
  df_sum_sp_year_mean <- as.data.frame(cbind(rep(species,length(yearRange)), rep(crop_type,length(yearRange)), yearRange, 
                                             matrix(out_sum_sp_year_mean, nrow = 30, byrow = TRUE)))
  colnames(df_sum_sp_year_mean) <- c("Crop","CropType","Year","Yield","Biomass","Duration")
  df_sum <- rbind(df_sum, df_sum_sp_year_mean)
}
df_sum <- df_sum[-1,]
df_sum$Year <- as.numeric(df_sum$Year)
df_sum$Yield <- as.numeric(df_sum$Yield)
df_sum$Yield[is.nan(df_sum$Yield)] <- NA
df_sum$Biomass <- as.numeric(df_sum$Biomass)
df_sum$Biomass[is.nan(df_sum$Biomass)] <- NA
df_sum$Duration <- as.numeric(df_sum$Duration)
df_sum$Duration[is.nan(df_sum$Duration)] <- NA
# converting a list-column to a regular vector column
df_sum$CropType <- factor(unlist(df_sum$CropType))
df_sum$Crop <- factor(unlist(df_sum$Crop))

p <- ggplot(df_sum, aes(Duration, Biomass, color = CropType)) +
  geom_point(aes(size = Yield, frame = Year, ids = Crop)) +
  scale_x_log10()

ggplotly(p)
saveWidget(p, file=paste0("./Simple_output_discover/26crops/country_figures/Trace_Animations.html"))